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ABSTRACT 


Based on a hybrid galactic cosmic-ray transport model, which incorporated MHD global heliospheric data into 
Parker’s cosmic-ray transport equation, we studied the behavior of the transport of galactic cosmic rays and the 
corresponding gradients in their flux near the heliopause (HP). We found that, (1) by increasing the ratio of the 
parallel diffusion coefficient to the perpendicular diffusion coefficient in the interstellar magnetic field of the outer 
heliosheath, the simulated radial flux near the HP increases as well. As the ratio multiplying factor reached 10'°, 

the radial flux experienced a sudden jump near the HP, similar to what Voyager 1 observed in 2012. (2) The effect 
of changing the diffusion coefficients’ ratio on the radial flux variation depends on the energy of the cosmic rays, 
the lower the energy, the more pronounced the effect is. (3) The magnitude of the diffusion coefficients also affect 
the radial flux near the HP, the modulation beyond the HP varies by adjusting the magnitude multiplying factor. 
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1. INTRODUCTION 


After nearly four decades since it was lunched, Voyager 1 is 
now more than 130 AU from the Earth. Recent observations 
indicate that Voyager 1 may have already entered into the local 
interstellar medium (ISM). It was found that the above 70 MeV 
galactic cosmic-ray intensity increased about 30% on 2012 
August 25 as the spacecraft was at 121.7 AU, and at the same 
time, the anomalous cosmic-ray intensity detected by the Low 
Energy Charged Particle instrument decreased by an order of 
magnitude (Krimigis et al. 2013; Stone et al. 2013; Webber & 
McDonald 2013). In addition to the cosmic-ray intensity 
change, Decker et al. (2012) found that the plasma speed at 
Voyager I is nearly zero after 2010 April. Since it is widely 
accepted that the galactic cosmic-ray intensity should increase 
while anomalous cosmic-ray intensity should decrease when 
crossing the heliopause(HP), these signatures are consistent 
with the HP crossing by Voyager 1. 

Burlaga et al. (2013), however, found that the magnetic field 
direction did not change significantly during these cosmic-ray 
intensity changing events, and that the magnetic field lines still 
coincided more or less with the overall structure of the 
heliospheric spiral line. This is far different from the direction 
of the expected ISM magnetic field even after a draping of field 
lines by the center heliosheath has been taken into account. At 
first, this caused some uncertainty; thus, it was suggested that 
Voyager I had crossed a well-defined boundary for energetic 
particles that was possibly related to the HP (Webber & 
McDonald 2013), and that Voyager J was in a “heliosheath 
depletion region” (Burlaga et al. 2013; Krimigis et al. 2013; 
Stone et al. 2013). More recently, the plasma wave instrument 
on Voyager I (Gurnett et al. 2013) detected locally generated 
plasma oscillations with a frequency consistent with the plasma 
density of the local ISM. Thus, it has become more certain that 
Voyager I is in the local ISM, or at least in the very local ISM. 
Despite the latter observation, some debate still continues, e.g., 
Fisk & Gloeckler (2014) proposed a model that is consistent 


Sun: heliosphere 


with all of the Voyager J observations, but assuming that 
Voyager 1 is still in the inner heliosheath. Later, Gloeckler & 
Fisk (2014) even provided a test for this model. Meanwhile, 
Borovikov & Pogorelov (2014) argued that Voyager 1 might 
have been inside eddies formed by plasma instabilities at the 
HP. Alternative arguments were presented by e.g., Grygorczuk 
et al. (2014) and Strumik et al. (2014). Future data from 
Voyager I and the expected crossing of the HP by Voyager 2, 
perhaps sooner than expected, will surely enlighten us. Burlaga 
& Ness (2014) showed that the sector boundary predicted 
by Fisk & Gloeckler (2014) had not been observed through 
2014/151 and that the field direction, but not the magnitude, 
has been quite constant for this period. In addition, both the 
larger than 70 and 0.5 MeV cosmic-ray counting rates detected 
by Voyager I have been essentially constant for nearly three 
years (http://voyager.gsfc.nasa.gov/heliopause/data.html). 

The observational data from Voyager 1 have stimulated also 
several theoretical investigations of galactic cosmic-ray trans- 
port near the HP. Scherer et al. (2011) and Strauss et al. (2013) 
argued that the HP is not the modulation boundary for galactic 
cosmic rays so that there should be some level of modulation 
that happened beyond the HP (in the outer heliosheath (OHS)). 
On the other hand, Kota & Jokipii (2014) arrived at the opinion 
that galactic cosmic ray modulation is very small beyond the 
HP if the diffusion coefficients in this region are set to be large 
enough. Later, Guo & Florinski (2014) shared the same 
opinion that galactic cosmic-ray modulation beyond the HP is 
negligible. 

In such an atmosphere where different opinions exist, we are 
motivated to perform an independent study on this issue and 
strive to contribute some understanding for the community. 
Specifically, we will use a stochastic differential equation 
(SDE) numerical approach to investigate the galactic cosmic- 
ray transport in a global heliosphere from an MHD-neutral 
kinetic simulation, which is thought to be as realistic as 
possible. A brief description of our numerical model, the hybrid 
transport model, will be presented in the subsequent section. In 
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Figure 1. MHD simulated profile of the magnetic field magnitude along the 
meridional plane (X-Z plane). Trajectories of the two spacecraft Voyager 1 and 


Voyager 2 are projected onto the same plane as shown by the white lines. The 
black curves outline the profiles of the termination shock and HP, respectively. 


the third section, we will discuss our simulation results: the 
galactic cosmic-ray spectra, the simulated radial flux variations 
by modifying the diffusion coefficients differently, and some 
possible mechanism for the simulation results. We will 
conclude in the last section. 


2. NUMERICAL METHOD 


We performed this investigation by using a hybrid galactic 
cosmic-ray transport model, incorporating the output from a 
global heliospheric MHD model into the galactic cosmic-ray 
transport SDE-type code. See Luo et al. (2013) for details of 
this approach. 


2.1. Realistic MHD Global Heliosphere Model 


The numerical MHD global heliospheric data, which supply 
the plasma and magnetic field background to the galactic 
cosmic-ray transport model, is obtained by solving a set of 
MHD equations which describe the interaction between the 
ISM and the solar wind plasma flow. The ISM is partially 
ionized. Neutral ISM atoms interact with plasma through 
charge exchange and photoionization provide an extra source 
of particle momentum and energy. A multi-component model 
of neutral atoms is used where the latter are subdivided into 
populations depending on the place of their origin and further 
treated gasdynamically as separate fluids (Zank 1996; Pogor- 
elov et al. 2006). For details, see Pogorelov et al. (2009a, 
2009b, 2012, 2013). MHD simulations are performed on a 
grid, while their results may be required at arbitrary points 
inside the computational regions. These are obtained by 
interpolation. 

Figure | shows a snapshot distribution of the magnetic field 
in the meridional plane formed by the Sun’s rotation axis (the 
Z-axis) and the unperturbed LISM velocity vector. As a result, 
the X-axis belongs to the meridional plane and is directed 
upwind from the LISM. Similarly, Figure 2 shows the snapshot 
distribution of the plasma speed’s radial component in the 
meridional plane. The trajectories of the two Voyager space- 
craft, the profiles of the termination shock and HP are also 
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Figure 2. MHD simulated solar wind radial component along the X-Z plane. 


shown in Figure |. Based on these two figures, we note that 
there is a strong V/—-V2 asymmetry of the heliopause, which is 
created by the interstellar magnetic field (Pogorelov et al. 
2008). It should be stated that within the confines of our MHD- 
neutral model, the position of the heliopause in the VZ 
trajectory direction is overestimated (146 AU instead of the 
observed 122 AU). However, this does not cause many 
qualitative differences. 


2.2. Hybrid Galactic Cosmic-Ray Transport Model 


Our investigation is based on the Parker’s transport equation 
(Parker 1965), which contains the solar wind velocity V, the 
averaged drift velocity (Vp), a diffusion term V - (K® - Vf) 


and adiabatic energy changes ECV - V)Əf/ð In p, with p as the 
momentum and f as the cosmic-ray distribution function. 


vv + (Vo)) - Vf 


1 of 
.(K®. = ; : 1 
+V- (K®. Vf) + S Yous (1) 


In the magnetic field coordinate system, the diffusion tensor 
K can be written as 


Ky O O 
K®=10 kı O|. (2) 
0 0 KL 


Following previous work (Luo et al. 2011, 2013), we adopt 
the following “traditional” forms for these diffusion coeffi- 
cients: 


0.5 
Ky = nobl 2 pe Ga) 
0.5 
Be 
61 = Kio8|+| =, (3b) 
0 B 


where, & 10 and kj are both constants, B is the magnetic field 
magnitude, and ( is the ratio of particle speed to the speed of 
light. The po parameter is a reference momentum (in our case 
1 GeV/c) and Beq is the magnitude of the heliospheric magnetic 
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field at the Earth (in the heliospheric equator at 1 AU). Using 
the MHD plasma data, we can calculate the diffusion tensor at 
any point X for particle with a momentum of p. As for the 
average drift velocity, we use the classical form (Jokipii & 
Thomas 1981; Potgieter & Moraal 1985) 


mafri, (4) 
3q B? 
As for the diffusion coefficients, it depends on the spatial 
dependence of the magnitude of the heliospheric magnetic 
field, but also on the direction of this field. To obtain the 
convection and adiabatic energy change terms is straightfor- 
ward, since they are related to the plasma speed profile, the 
latter on the divergence of this profile. 
Following Markov’s stochastic method (Zhang 1999), the 
transport equation can be changed to the following time- 
backward SDEs. 


dX = (V - K® — V — (Vp))ds + Y aodW (s), Sa) 


dp = Zp - V)ds. (5b) 


In this equation, dW, (s) is the Wiener process, and it can be 
generated in each step using a Gaussian distribution random 
number. Based on this method, (X, p) constitutes the phase 
space for the distribution function f. We set |X| < 300 AU, the 
polar angle 0 € [0, 7), the azimuthal angle y € [0, 27), and 
momentum p € (0, 1000p); po is the initial momentum for 
tracing. In order to get the value of f (Xo, po) at the point 
(Xo, Po) in phase space, we ran a large number of stochastic 
trajectories from (Xo, pọ) backward in time until they hit the 
modulation boundary for the first time. A similar approach was 
followed by Kopp et al. (2012). 

The solution for the modulated cosmic-ray distribution 
function is 


SE. p) = (fo (Xe r)) (6) 


where fẹ is the boundary condition where the stochastic 
trajectories hit the boundary for the first time; ( } denotes 
the ensemble average. Each stochastic trajectory represents 
a number of pseudo-particles proportional to the boundary 
value. 


3. RESULTS AND DISCUSSION 


In the following, we present our simulation results using this 
hybrid model. Specifically, we first simulate the galactic 
cosmic-ray spectrum at different radial locations to test and 
validate our code, then we modify the diffusion coefficients 
beyond the HP to explore how changing their ratio and their 
absolute values could affect the cosmic-ray transport ahead and 
beyond the HP. 


3.1. Cosmic-Ray Spectra 


We first run a series of test simulations for cosmic-ray 
spectra along the Voyager 1 direction at different locations. The 
interstellar spectrum is specified at 300 AU, which is our 
simulation boundary. We take note that, recently, the very local 
interstellar spectrum (LIS) for protons has been newly 
determined (Webber et al. 2013; Potgieter 2014), but since 
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Figure 3. Simulated galactic cosmic-ray spectra along Voyager 1’s trajectory at 
different radial locations. The unit for flux is arbitrary for Figures 4, 6, and 7. 


the radial gradient of the flux is not affected by the exact LIS 
spectral shape (Luo et al. 2013), we still adopt the previously 
used form in this study: 


fam P) x (mê + p?) p, (7) 


where mp is the mass of the proton. 

Since our study is mainly about the galactic cosmic-ray 
radial gradient, its absolute level is not crucially important. The 
simulation results shown below are therefore in relative 
intensity with arbitrary units of the flux j. 

Figure 3 illustrates these simulation results, which clearly 
demonstrate galactic cosmic-ray modulation and its basic 
features from the modulation boundary to the inner helio- 
sphere. The simulations were computed by setting 
Kjo = 50 x 10% cm? s™' and x1, = 5 x 10% cm? s™' in Equa- 
tion (3), so that the ratio of the diffusion coefficients is 
K/K = 10 in the whole simulation domain. It should be 
noted that the spectrum at 155 AU is lower than the interstellar 
spectrum, because this particular choice of modulation 
parameters causes modulation in this region, which is beyond 
the HP. This issue is further addressed below since it is 
presently debated as mentioned above. Based on current 
understanding of the ISM (Armstrong et al. 1995; Büsching & 
Potgieter 2008; Shalchi et al. 2010), the diffusion there is quite 
different from the situation inside the heliosphere. In the 
following, we will investigate how the variation of the 
diffusion beyond the HP affects the cosmic-ray transport there 
by using our numerical approach. 


3.2. The Radial Intensity Variations 


As suggested by Büsching & Potgieter (2008) and Shalchi 
et al. (2010), cosmic-ray propagation in the ISM is quite 
different because of the properties of magnetic turbulence 
inside the ISM. The perpendicular diffusion coefficient is 
assumed to be much smaller than the parallel diffusion 
coefficient in the ISM. Calculations based on interstellar 
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Figure 4. Simulated radial flux for 100 MeV protons. The three curves illustrate three different cases by changing the ratio of the parallel diffusion coefficient to the 
perpendicular diffusion coefficient in the OHS: (A) blue solid curve with 4 / « = 10!°, (B) dashed brown curve with ky / Kı = 104, and (C) dotted purple curve with 
Kj /K 1 = 10. In the rest of the heliosphere, Rij /£ 1 = 10 for all cases, which means that for the “no change” scenario, this ratio remains the same. 


diffusion models by, e.g., Strong et al. (2007) indicate that the 
scale of the diffusion coefficient is on the order of 1078 cm? s~. 
This level of diffusion is mostly from particle diffusion along 
the local magnetic field direction, or «j. As for the OHS, the 
situation is still unclear, but we anticipate that parallel diffusion 
is also very effective in the OHS. The magnetic field turbulence 
in the OHS should be quite small as measured by Voyager 1 
(Burlaga et al. 2014) and inferred from the IBEX ribbon 
(Gamayunov et al. 2010). A quiet magnetic field warrants large 
parallel diffusion and small perpendicular diffusion, so we 
expect k/k to be significantly large in the OHS. 

We investigate the diffusion coefficients in the OHS by 
adjusting the ratio of the parallel diffusion coefficient to the 
perpendicular diffusion, as was done by Strauss et al. (2013), 
as well as their absolute values. In this numerical approach, we 
also want to illustrate and understand how cosmic-ray 
modulation behaves inside of the HP (upstream), closer and 
across the HP when these type of changes are made. 

Figure 4 shows the simulated radial flux for 100 MeV 
protons along Voyager ls direction (polar angle 0 = 56°, 
azimuthal angle ¢ = 4°). The three curves represent three 
different cases: (A) for the blue solid curve the ratio of Ky /K 
was increased by 10'° in the OHS. (B) For the brown dashed 
curve, the ratio of kj /x, was magnified by a factor of 104 in the 
OHS. (C) For the purple dotted curve, we keep the ratio as a 
constant in the simulation domain. For these scenarios, 
excluding the OHS, we set the ratio of «j/x = 10 everywhere 
inside the heliosphere, that is, in all upstream regions from the 
HP. The details of the value of the diffusion coefficients and 
magnetic field magnitude variations along the Voyager 1 
direction in the simulation domain are shown in Figure 5. In the 
outer heliosphere, ,, is very close to x]. It increases beyond 
the HP as we set the x| /%_ ratio to 10!° and it reaches the value 
of « after crossing the HP. The magnetic field magnitude 
decreases inside the supersonic solar wind region, such as in 


Parker’s interplanetary magnetic field model. It increases after 
crossing the termination shock, probably due to the shock 
compression. Around 130 AU, because of the current sheet 
crossing, it decreases again, and a “magnetic wall” with 
magnitude increase can be seen around 150—180 AU. 

From Figure 4, it follows that the corresponding radial flux 
gradient becomes significantly different after adjusting the ratio 
of the parallel diffusion coefficient to the perpendicular 
diffusion coefficient. The higher the ratio, the larger the radial 
gradient near the HP. As the ratio approaches 10!”, like the K, 
curve trend, the radial gradient reaches very large values as the 
flux jumps to the interstellar value in a very short distance; 
which is quiet similar to what Voyager J observed in 2012 
August (Webber & McDonald 2013). In addition, beyond the 
HP, the simulations demonstrate that the flux and correspond- 
ing gradient differ significantly depending on the assumed 
ratios. If the ratio is unchanged, the modulation simply 
continues in the OHS as if it is part of the global heliospheric 
medium. 

We also expanded our simulation for protons with different 
energy. Figure 6 shows the simulated radial flux for 200, 100, 
80, and 50MeV, which is above the anomalous proton 
energies. We set j/k] = 10'° in the OHS for all of these 
simulations. Similar to case (C) of Figure 4, the radial flux 
jumps upward near the HP for all of these energies. For 
50 MeV, this jump near the HP contributes about 25% of total 
modulation; for 80 MeV, the jump contributes 17% of the total 
modulation; for 100 MeV, it contributes 15% and for 200 MeV 
protons, only 12%. Evidently, as the energy increases, the 
radial flux jump level near the HP became less. As a result, the 
effect of the ratio variation becomes less and the increasing 
factor becomes less important near the HP for higher energy 
protons, but clearly quite significant at lower energies. 

We also perform simulations by increasing the value of each 
individual diffusion coefficient inside the OHS while keeping 
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Figure 5. Values of the different diffusion coefficients and magnetic field magnitude, as indicated in the legend, along the Voyager 1 direction as a function of radial 


distance in the simulation domain for cases (A) and (C). 
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Figure 6. Simulated radial flux of galactic protons as a function of radial distance for four different proton energies. The emphasis here is on how the flux changes with 


decreasing energy across the HP. 


the ratio the same as used inside of the HP (upstream). In 
Figure 7, the results are shown together with two scenarios 
from Figure 4, repeated as references (the black dotted and 
dashed curves). We multiply both parallel and perpendicular 
diffusion coefficients by a factor of five in the OHS. The flux 
(solid blue curve) increases correspondingly, but the value 
around the HP is still lower than the interstellar value, with 
some modulation still occurring beyond the HP for this case. 
As we set the multiplying factor for both diffusion coefficients 
to 100, the flux (dash-dotted curve) value around the HP 
reaches the interstellar value very quickly, with the jump in the 
flux at the HP not as obvious as before. This simulation shows 


that either changing the ratio or the value of the individual 
diffusion coefficients affects galactic cosmic-ray transport near 
the HP. It appears that the observed jump in the flux at the HP 
and a case of no modulation beyond the HP, requires a large Ky 
and fj /K, ratio. 

It is worth mentioning here that, based on the stochastic 
method as utilized here, we are able to trace individual pseudo- 
particle trajectories in the simulation domain. Because the 
pseudo-particles have the same distribution as real particles 
entering at the modulation boundary, for a case of little 
modulation near the HP region, we can approximately consider 
these pseudo-particles as real particles. Since this requires a 
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Figure 7. Simulated radial flux for 100 MeV protons across the HP. Instead of changing the ratio of the diffusion coefficients, the individual diffusion coefficients are 
now changed as indicated in the legend while keeping the ratio the same as used inside (upstream) of the HP. 


lengthy description, we refrain from showing such trajectories 
in this manuscript but as a next step, we plan to investigate the 
exiting characteristics of these pseudo-particles, which is where 
real particles are entering the heliosphere in order to reach the 
Voyager 1 position. Hopefully, this will yield some under- 
standing of the possible physical processes throughout the 
heliosphere, also beyond the HP, and from where cosmic-ray 
particles can actually be transported before reaching Voyager 1 
and Voyager 2. 


4. SUMMARY 


In this paper, by incorporating the output of a global MHD 
heliospheric model into the galactic cosmic-ray transport 
model, we constructed a hybrid cosmic-ray transport model. 
Based on this model, we investigated the behavior and features 
for cosmic-ray transport near the HP. We presented proton 
fluxes showing that the radial flux near the HP can already be 
modulated by the OHS if the diffusion coefficient ratio Ky /K is 
set to a small value. By adjusting this ratio to a very large value 
in the OHS, it was found that radial flux exhibits a sudden 
upward jump near the HP, which is similar to what Voyager 1 
observed in 2012. Similar features have also been shown by 
Strauss et al. (2013), Guo & Florinski (2014), and Kota & 
Jokipii (2014). Modulation beyond the HP seems indeed 
possible, but since we do not know the exact values for the 
relevant diffusion coefficients it is difficult to predict how large 
this modulation may be. We found that the effect of changing 
the ratio on the jump in flux is closely related to the energy of 
the protons, the lower the energy, the larger the effect. After 
adjusting the magnitude of the individual diffusion coefficients, 
the radial flux also differs significantly. However, this does not 
give a significant jump of flux at the HP. 

We also showed that there is little modulation occurring 
beyond the HP after multiplying the values of the individual 
diffusion coefficients by a small factor, while K/K remains 
the same as it is inside of the HP, upstream toward the Sun. 


However, at this stage, without published observational 
cosmic-ray data, it is difficult to figure out if this scenario is 
plausible. For future work, we plan to investigate this further 
and to link the Voyager 1 observations with a realistic physical 
environment near the HP, thus constraining the range of 
relevant parameters. This should contribute to a further 
understanding of the recent observations by Voyager I and 
what it may imply for Voyager 2. In this study, we used the 
analytic forms as done before (Zhang 1999) for the 
perpendicular and parallel diffusion coefficients, in particular, 
using a simple rigidity dependence, which, if changed, could 
affect the results in terms of energy dependence shown here. In 
a next paper, changing this to more complex forms will be 
investigated. 
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